Microscopic universality with dynamical fermions 
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It has recently been demonstrated in quenched lattice simulations that the distribution of the low- 
lying eigenvalues of the QCD Dirac operator is universal and described by random-matrix theory. 
We present first evidence that this universality continues to hold in the presence of dynamical 
quarks. Data from a lattice simulation with gauge group SU(2) and dynamical staggered fermions 
are compared to the predictions of the chiral symplectic ensemble of random-matrix theory with 
massive dynamical quarks. Good agreement is found in this exploratory study. We also discuss 
implications of our results. 
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It was conjectured a few years ago that the mi- 
croscopic spectral properties of the QCD Dirac opera- 
tor, in particular the distribution of the low-lying eigen- 
values, are universal and can be obtained in effective 
theories that are much simpler than QCD. This con- 
jecture rests on the observation j^] that in the range 
1/A <IC I/i/'* ^ l/»Ti7r, where A is a typical hadronic 
scale, V is the space-time volume, and is the pion 
mass, the mass dependence of the finite- volume partition 
function of QCD is completely determined by global sym- 
metries. Thus, the spectral properties of the Dirac oper- 
ator on the "microscopic" scale ~ where E is the 
absolute value of the chiral condensate, can be computed 
in effective theories with only the global symmetries as 
input, such as the framework of effective Lagrangians 0] 
or chiral random- matrix theory (RMT) ||]. 

The global spectral density of the Euclidean Dirac op- 
erator = + ig4 is given by p(A) = {Y^n^i^ " A„))a, 
where the A„ are the eigenvalues of «0 and the aver- 
age is over all gauge field configurations A weighted by 
exp(— S'qcd)- The Banks-Casher formula, E = TTp{0)/V 
[Q, relates the spectral density at zero virtuality to the 
chiral condensate. If chiral symmetry is spontaneously 
broken, this relation implies that the spacing of the low- 
lying eigenvalues is ~ 1/{V'E), rather than ~ as 
in the case of the non-interacting Dirac operator. Thus, 
the distribution of the low-lying eigenvalues is of great 
interest for a better understanding of the phenomenon of 
spontaneous chiral symmetry breaking. 

It has recently been shown that the microscopic spec- 
tral properties of the staggered lattice Dirac operator in 
quenched SU(2) are indeed universal and described by 
chiral RMT [£|j6|. However, for a better understanding 
of hadronic properties it is important to go beyond the 
quenched approximation. One of the problems that arises 
in unquenched lattice simulations is that one cannot go 
to arbitrarily small quark masses. An important point 
is how small the masses of the dynamical quarks have to 
be so that they are really dynamical, i.e., lead to results 
that arc different from those obtained in the quenched 
approximation. In this context, the main question that 



will be addressed in this work is the following. What is 
the effect of light dynamical quarks with masses of order 
^ 1/{VY,) on the low-lying spectrum of the Dirac opera- 
tor, and can this effect be described by results obtained 
in chiral RMT? The answer to these questions together 
with the analytical information available from RMT leads 
to a better understanding of the mass and energy scales 
in the problem and has important technical implications 
with regard to extrapolations to various limits (chiral, 
thermodynamic, continuum) that are difficult to take on 
the lattice. Moreover, analytical knowledge of the dis- 
tribution of the smallest eigenvalues may be interesting 
from an algorithmic point of view. 

We will mainly be concerned with the distribution of 
the smallest positive eigenvalue, P(Amin), and with the 
microscopic spectral density, Ps{z), defined by |^ 



Ps{z) 



1 



lim n I 



(1) 



This definition amounts to a magnification of the region 
of low-lying eigenvalues by a factor of VE and leads to 
the resolution of individual eigenvalues. The claim is 
that the quantities P{Xadn) and Ps{z) are universal (see 
Ref. 1^ for a summary of existing evidence) and can be 
computed exactly in an effective theory. We will con- 
centrate on chiral RMT although identical results could 
also be obtained from the finite- volume partition function 
computed from an effective Lagrangian 

Let us recall how the random-matrix model is con- 
structed. In a chiral basis, the weight function of full 
QCD with Nf fiavors can be written as 
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where Sg\ is the gluonic action, T is a matrix representing 
i0, and the m/ are the masses of the dynamical quarks. 
We now replace T by a random matrix W with N rows 
and N + v columns (|j^| <C N). The dimensionless space- 
time volume can be identified with 2N , and v plays the 
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role of the topological charge. The symmetry properties 
of W depend on the gauge group and on the represen- 
tation of the fermions |9|. We will be concerned with 
staggered fermions in SU(2) for which the elements of W 
are real quaternions. This is the chiral symplectic ensem- 
ble (chSE). The average over gauge field configurations 
is replaced by an average over random matrices, i.e., the 
factor e""^*' is replaced by a convenient distribution of 
the matrix W, often a simple Gaussian, but the results 
are insensitive to this choice PQ-|lJ]. For the present 
work, the second factor in Eq. is more interesting. In 
terms of the eigenvalues of W, it reads for the chSE M 
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in a sector of topological charge i/, where A is the Vander- 
monde determinant. Thus, in the random-matrix model 
the fermion determinants can be taken into account with- 
out further assumptions. From Eq. it is intuitively 
clear that the presence of the fermion determinants will 
affect the microscopic spectral quantities only if the m f 
are on the order of the smallest eigenvalue or, in other 
words, on the order of the mean level spacing near zero 
Thus, we require m/ ^ to observe an ef- 

fect on P(Ainin) and Ps{z)- For quark masses much larger 
than this, we should simply obtain agreement with the 
RMT-results computed in the quenched approximation. 
This was already observed in Ref. where the lattice 
data of Ref. were analyzed. 

RMT-predictions for the microscopic spectral quanti- 
ties in the presence of massive dynamical quarks have 
recently been computed for the chiral unitary ensem- 
ble (see also 1 13 ) and the unitary ensemble . 
However, we want to compare lattice data with random- 
matrix predictions for the chSE. In this case, closed ana- 
lytical expressions are at present only known in the chiral 
limit p9|-pl|. Analytical work for the massive case is in 
progress, but in the meantime we have computed the 
RMT-predictions for the chSE with massive quarks nu- 
merically using the methods of Ref. |Q . Essentially, one 
has to do an iterative computation of skew-orthogonal 
polynomials which obey orthogonality relations deter- 
mined by a weight function involving the fermion deter- 
minants. To avoid cancellation problems, we have used 
a multi-precision package 1 22 1 . The calculation was done 
with matrices W of finite dimension up to A'' = 100. 
The remaining finite- A^ effects can be estimated from a 
calculation in the chiral limit and are on the order of 1%. 

We now summarize some RMT-results for the chSE in 
the chiral limit which we will need in the following. With 
a = Nf + 2\vl we have [|o|j23|J|] 



Ps{z) = z [Jl{2z) ~ J„+i(2z)J„_i(2z)] 
1 f^' 
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where J denotes the Bessel function. Rescaling A„ 
VT, as in Eq. (0), we have for Nf = ly = ^ 
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where / is the modified Bessel function. Further analyti- 
cal results for P(Ainin) are known if a is odd, i.e., for odd 
Nf. For example, we have for Nf — 1 and v — pl| 



-P(Amin) — 2'^™'"P^2(2Aniin) ^ ^o(2Annn) + 1] 6 



(6) 



We shall also need P{Xmin) for Nf = 2 and 4 which in the 
absence of closed analytical expressions can perhaps most 
easily be obtained from results computed by Kaneko . 
The essential ingredients are the zonal polynomials at the 
symmetric point. We obtain after some algebra 
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where T{x) = 1 + X^tfci '^d x'^ with 
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Here, k denotes a partition of the integer d, 1{k) its 
length, |k| its weight, and k' the conjugate partition. 
In Eq. (@), a partition k is identified with its diagram, 
K = {s = 1 < i < /(k), 1 < j < Ki}. The Taylor 

series for T{x) is rapidly convergent, and the curves for 
Nf = 2 and 4 can easily be computed to any desired 
accuracy. 

The RMT-results for P(Ainin) and Ps{z) are sensitive 
to the topological charge Thus, they are only univer- 
sal in a given sector of topological charge. In the con- 
tinuum limit, the overall result for a given quantity is a 
weighted average over these sectors. RMT can predict 
the quantity for definite v but cannot predict the corre- 
sponding weights. In Ref. ||] it was found that for values 
of /3 = 4/(7^ up to 2.4 the quenched SU(2) lattice data 
were consistent with v = 0. Here, we only consider val- 
ues of f3 in the strong-coupling region, thus everything 
should be described by the RMT-results for i/ ~ 0. We 
will not address the question of topology in the following. 

Our lattice simulations were performed with gauge 
group SU(2) and staggered fermions on an S'* lattice. The 
boundary conditions were periodic for the gauge fields 
and periodic in space and anti-periodic in Euclidean time 
for the fermions. A hybrid Monte Carlo algorithm ^5 26 
was used to generate a large number of independent con- 
figurations. For the diagonalization of —0^ we have 
employed the CuUum-Willoughby version of the Lanc- 
zos algorithm |M. The complete spectra where checked 
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TABLE I. Simulation parameters and observables of our 
lattice calculations on an S** lattice. {Tpip) (m) and E are given 
in units of a~'^. 






Nf 


ma 


conf. 


{ipip)im) 


S 




1.8 


4 


0.015 


1799 


0.1517(14) 


0.1250(46) 


7.68 


1.3 


8 


0.0075 


1023 


0.2669(12) 


0.2846(38) 


8.74 


1.3 


8 


0.0055 


1042 


0.2470(11) 


0.2800(45) 


6.31 


1.3 


8 


0.00375 


703 


0.2073(22) 


0.2772(42) 


4.26 



against the analytical sum rule tr(— 0^) = V for the dis- 
tinct eigenvalues of |^^. All physical quantities will 
be quoted in units of the lattice spacing a, in contrast to 
Ref. where the unit was 2a. 

The RMT-results contain one single parameter VH 
which is determined by the lattice data via the Banks- 
Casher relation, VH = 7rp(0). (Note that p{X) is nor- 
malized to V .) Thus, the random-matrix predictions are 
parameter-free. The simulation parameters and observ- 
ables are summarized in Table |. We also give values for 
{^pip){m) = {m/V)J2x„>oi^l 4- m2)-i. Note that we 
have /i — mVYi ^ 0{\) as desired, where by 0(1) we 
mean much smaller than 8'^ = 4096. 

Our results for (3 = 1.8 and Nf = 4 flavors of mass 
ma — 0.015 are plotted in Fig. |l]. Note that we have to 
compare the data with the RMT-predictions for N f — 2 
in Eq. (||) for the following two reasons. First, because 
of the finite lattice constant the original U(4) (g) U(4) chi- 
ral symmetry of the (local) lattice action is broken down 
to 15(1) (8) U(l). This would imply to use Nf = I in 
Eq. (||). Second, because of a global charge conjuga- 
tion symmetry which is special to SU(2) all eigenvalues 
are twofold degenerate. Equation (^) was derived for 
fermions in the adjoint representation, and the eigenval- 
ues in the fermion determinants of Eq. (^ are assumed 
to be non-degenerate |9|. Therefore, the number of fla- 
vors must be doubled. Note that this is only an issue for 
SU(2). For these two reasons, we necessarily have to use 
Nf = 2 in the random-matrix results. 

After these comments, let us discuss Fig. ||. The pres- 
ence of the fermion determinants leads to a suppression 
of small eigenvalues. The fact that the lattice data dif- 
fer from the RMT-curves for the quenched approxima- 
tion [Nf = 0, or /i — oo) means that the dynamical 
quarks were light enough to affect the microscopic spec- 
tral quantities. However, the data are still closer to the 
RMT-curves for Nf — than to the ones for the chiral 
limit (^ = 0). This means that the dynamical quarks 
with a rescaled mass of /i = 7.68 are still quite heavy as 
far as the microscopic spectral properties are concerned. 
We will explore smaller masses in Fig. I below. The lat- 
tice data agree reasonably well with the random-matrix 
predictions for Nf = 2 with the appropriate value of /i. 
The systematic deviations and the question of the contin- 
uum limit will be discussed below and in the conclusions. 
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FIG. 1. Distribution of the smallest eigenvalue (left) and 
microscopic spectral density (right) of the lattice Dirac oper- 
ator for the simulation parameters indicated above the figures. 
The histograms represent the lattice data. The full lines are 
the RMT-predictions for Nf ^ 2 and /i = mVT, = 7.68, the 
dashed lines those for the quenched approximation, and the 
dotted lines those for Nf = 2 in the chiral limit, respectively. 
Note that there is no free parameter involved. 

For a larger number of flavors, the curves for the 
quenched approximation and those for the chiral limit are 
farther apart from each other so that the influence of mas- 
sive dynamical quarks can be seen more easily. Therefore, 
we have also performed lattice simulations with Nf = 8 
flavors, shown in Fig. ||. A value of /3 = 1.3 was chosen 
to ensure that we are in the broken phase for all values 
of the bare quark mass |^6| . As expected from the above 
discussion, we have to compare the lattice data with the 
RMT-results for Nf = 4 and the appropriate value of 
mVTj. We see in Fig. || that, as the mass of the dynam- 
ical quarks is lowered, the lattice data move away from 
the quenched curves towards the RMT-curves computed 
in the chiral limit. This is precisely what we expect. Let 
us discuss the numbers. Our qualitative criterion was 
^ - 0(1). For ^ = 8.74 (the top two figures in Fig. |), 
we are still rather close to the quenched RMT-curves. 
For /i — 4.26 (the bottom two figures in Fig. H), we are 
about half-way in between the quenched approximation 
and the chiral limit. We can now say more quantitatively 
that values of fj, around 10 yield results which are close to 
those obtained in the quenched approximation, whereas 
values of smaller than 4 yield results which are closer 
to those obtained in the chiral limit. The mass scale on 
which dynamical quarks affect the microscopic spectrum 
of the Dirac operator is thus identified quantitatively. 

In Fig. ^, we compare the data for P(A,„in) computed 
for Nf = 8 and the lightest quark mass, m = 0.00375, 
with the RMT-predictions for various numbers of flavors 
and the appropriate value of /i. The fact that the data 
agree well with the RMT-result for Nf — 4 supports the 
two arguments made in connection with Fig. ^ In the 
continuum limit, one would expect agreement with the 
RMT-curve for Nf = 16. 

The agreement between the various lattice data and 
the RMT-predictions is not perfect. There are some sys- 
tematic deviations which we believe to be largely due 
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FIG. 2. Same as Fig. |l| but for different simulation param- 
eters (indicated above the figures). 



to the finite size of the lattice. Similar finite-size effects 
were observed in Ref. |^] . The present exploratory study 
was restricted to a relatively small lattice size since a 
large number of independent configurations is needed. 
One also observes from the figures that the agreement 
with RMT is better for smaller values of /x, as expected 
from the discussion in the introduction. While it would 
be desirable to perform a comprehensive study for larger 
lattices and a number of different values of Nf, and m, 
the agreement between lattice data and random-matrix 
results seen in our figures is already very encouraging. 
We point out again that there is no free parameter in- 
volved. Another possibility would have been to fit the lat- 
tice data to the RMT-results by adjusting the parameter 
VT,. This is an alternative way of determining the chiral 
condensate which also seems to eliminate some finite-size 
effects In this case, the agreement in Figs. |l| and § 
would have been better. However, the primary purpose 
of this work is the demonstration of parameter-free agree- 
ment between lattice data and RMT-predictions. This is 
the reason why we determined E from p(0) and not from 
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FIG. 3. Distribution of the smallest eigenvalue for Nf = 8 
and m = 0.00375 (histogram) compared with the RMT-predi- 
ctions for various numbers of flavors (indicated in the figure) 
and /i = 4.26. 



a fit to RMT. 

To conclude, we comment on the implications of our 
finding that lattice data for the small eigenvalues of the 
Dirac operator are described by chiral RMT also in the 
presence of dynamical quarks of mass ~ 1/{VT,). First 
of all, as far as the microscopic spectral properties are 
concerned, the dynamical quarks are only really dynam- 
ical if their masses are not much larger than l/(yE). 
This provides information on the relevant mass scales 
in unqucnched lattice simulations. Once the applicabil- 
ity of chiral RMT to full QCD with massive dynamical 
quarks is firmly established, one can address practical 
applications. For quantities which are sensitive to the 
small eigenvalues, the analytical RMT-results can pro- 
vide guidance for extrapolations to the chiral limit. Also, 
as demonstrated in Ref. |2^, extrapolations to the ther- 
modynamic limit are facilitated since the RMT-results 
are derived in the limit V oo. Another interesting 
aspect which deserves further investigation is the con- 
tinuum limit where the original chiral symmetry of the 
action is restored and agreement with the RMT-results 
for the original number of flavors is expected. This might 
give an indication of how close to the continuum limit one 
actually is. However, much larger lattices are required 
to study such a transition. Moreover, we hope that the 
microscopic universality can perhaps be used in the de- 
velopment of hybrid fermionic algorithms that take into 
account the available analytical information. It would 
also be of great interest to extend our analysis to gauge 
group SU(3) where analytical RMT-results are already 
known. Finally, we hope to obtain analytical results for 
the chSE with massive dynamical quarks in the near fu- 
ture. 
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